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We describe here a library aimed at automating the solution of partial differential equations 
using the finite element method. By employing novel techniques for automated code generation, 
the library combines a high level of expressiveness with efficient computation. Finite element 
variational forms may be expressed in near mathematical notation, from which low-level code 
is automatically generated, compiled and seamlessly integrated with efficient implementations of 
computational meshes and high-performance linear algebra. Easy-to-use object-oriented interfaces 
to the library are provided in the form of a CH — h library and a Python module. This paper 
discusses the mathematical abstractions and methods used in the design of the library and its 
implementation. A number of examples are presented to demonstrate the use of the library in 
application code. 

Categories and Subject Descriptors: G.4 [Mathematical software]: Algorithm Design, Effi- 
ciency, User Interfaces; G.1.8 [Numerical analysis]: Partial differential equations — Finite Ele- 
ment Methods; D.f.2 [Programming techniques]: Automatic Programming 

General Terms: Algorithms, Design, Performance 
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1. INTRODUCTION 

Partial differential equations underpin many branches of science and their solution 
using computers is commonplace. Over time, the complexity and diversity of sci- 
entifically and industrially relevant differential equations has increased, which has 
placed new demands on the software used to solve them. Many specialized libraries 
have proved successful for a particular problem, but have lacked the flexibility to 
adapt to evolving demands. 

Software for the solution of partial differential equations is typically developed 
with a strong focus on performance, and it is a common conception that high 
performance may only be obtained by specialization. However, recent developments 
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Fig. 1. DOLFIN assembles any user-defined variational form on any (simplex) mesh for a wide 
range of finite elements using any user-defined or built-in linear algebra backend. DOLFIN relies 
on a form compiler for generation of the problem-specific code. The form compiler generates 
code conforming to the UFC (Unified Form-assembly Code) interface, either at compile-time or 
run-time, and the generated code is called during assembly by DOLFIN. A small interface layer 
is required for each linear algebra backend and is implemented as part of DOLFIN for PETSc, 
Trilinos/Epetra, uBLAS and MTL4. 



in finite element code generation have shown that this is only true in part [Kirby 
et al. 2005; Kirby et al. 2006; Kirby and Logg 2006; 2007]. Specialized code is still 
needed to achieve high performance, but the specialized code may be generated, 
thus relieving the programmer of time-consuming and error-prone tasks. 

We present in this paper the library DOLFIN which is aimed at the automated 
solution of partial differential equations using the finite element method. As will be 
elaborated, DOLFIN relies on a form compiler to generate the innermost loops of 
the finite element algorithm. This allows DOLFIN to implement a general and ef- 
ficient assembly algorithm. DOLFIN may assemble arbitrary rank tensors (scalars, 
vectors, matrices and higher-rank tensors) on simplex meshes in one, two and three 
space dimensions for a wide range of user-defined variational forms and for a wide 
range of finite elements. Furthermore, tensors may be assembled into any user- 
defined data structure, or any of the data structures implemented by one of the 
built-in linear algebra backends. For any combination of computational mesh, vari- 
ational form, finite element and linear algebra backend, the assembly is performed 
by the same code, as illustrated schematically in Figure 1, and code generation 
allows the assembly code to be efficient and compact. 

DOLFIN functions as the main programming interface and problem solving en- 
vironment of the FEniCS Project [FEniCS 2009], a collaborative effort towards 
the development of innovative concepts and tools for the automation of computa- 
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tional mathematical modeling, with an emphasis on partial differential equations. 
See Logg [2007] for a overview. All FEniCS components are released under the 
GNU General Public License or the GNU Lesser General Public License, and are 
made freely available at http://www.fenics.org. 

Initially, DOLFIN was a monolithic, stand-alone C++ library including imple- 
mentations of linear algebra, computational meshes, finite element basis functions, 
variational forms and finite element assembly. Since then, it has undergone a num- 
ber of design iterations and some functionality has now been "outsourced" to other 
FEniCS components and third-party software. The design encompasses coexistence 
with other libraries, and permits a user to select particular components (classes) 
rather than to commit to a rigid framework or an entire package. The design also 
allows DOLFIN to provide a complex and feature-rich system from a relatively 
small amount of code, which is made possible through automation and design so- 
phistication. For linear algebra functionality, third-party libraries are exploited, 
with a common programming interface to these backends implemented as part of 
DOLFIN. Finite element basis functions are evaluated by FIAT [Kirby 2009; 2004] 
and variational forms are handled by the Unified Form Language (UFL) library [Al- 
naes and Logg 2009; Alnaes 2009] and the FEniCS Form Compiler (FFC) [Logg et al. 
2009; Kirby and Logg 2006]. Alternatively, DOLFIN may use SyFi/SFC [Alnaes 
and Mardal 2009; 2010] for these tasks, or any other form compiler that conforms 
to the Unified Form-assembly Code (UFC) interface [Alnaes et al. 2009] for finite 
element code. Just-in-time compilation is handled by Instant [Alnaes et al. 2009]. 
FIAT, FFC, SyFi/SFC, UFC and Instant are all components of the FEniCS Project. 
Data structures and algorithms for computational meshes remain implemented as 
part of DOLFIN, as is the general assembly algorithm. 

Traditional object-oriented finite element libraries, including deal. II [Bangerth 
et al. 2007] and Diffpack [Langtangen 2003], provide basic tools such as compu- 
tational meshes, linear algebra interfaces and finite element basis functions. This 
greatly simplifies the implementation of finite element methods, but the user must 
typically implement the assembly algorithm (or at least part of it), which is time- 
consuming and error-prone. There exist today a number of projects that seek to 
create systems that, at least in part, automate the finite element method, including 
Sundance [Long et al. 2009], GctDP [Dular et al. 2009], FreeFEM++ [Pironneau 
et al. 2009] and Life [Prud'homme 2009; 2007]. All of these rely on some form of 
preprocessing (compilc-timc or run-time) to allow a level of mathematical expres- 
siveness to be combined with efficient run-time assembly of linear systems. DOLFIN 
differs from these project in that it relies more explicitly on code generation, which 
allows the assembly algorithms to be decoupled from the implementation of varia- 
tional forms and finite elements. As a result, DOLFIN supports a wider range of 
finite elements than any of the above-mentioned libraries since it may assemble any 
finite element variational form on any finite element space supported by the form 
compiler and finite element backend. 

The remainder of this paper is organized as follows. We first present a background 
to automated finite element computing in Section 2. We then present some general 
design considerations in Section 3 before discussing in more detail the design and 
implementation of DOLFIN in Section 4. We present in Section 5 a number of 
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Python code 

from dolfin import * 

mesh = UnitSquare (32 , 32) 

V = FunctionSpace(mesh, "CG", 1) 

v = TestFunction(V) 
u = TrialFunction(V) 

f = Expression("sin(x[0] )*cos(x[l] ) ") 

A = assemble (dot (grad(v) , grad(u))*dx + v*u*dx) 
b = assemble (v*f*dx) 

u_h = Functlon(V) 
solve(A, u_h. vector () , b) 

plot(u_h) 



Table I. A complete program for solving the reaction-diffusion problem (1) and plotting the 
solution. This and other presented code examples are written for DOLFIN version 0.9.5 (released 
in December 2009). 

examples to illustrate the use of DOLFIN in application code, which is followed by 
concluding remarks in Section 6. 

2. AUTOMATED FINITE ELEMENT COMPUTING 

DOLFIN automates the assembly of linear and nonlinear systems arising from the 
finite element discretization of partial differential equations expressed in variational 
form. To illustrate this, consider the reaction-diffusion equation 

-Au + u = f (1) 

on the unit square = (0, 1) x (0, 1) with f(x, y) = sin(x) cos(y) and homogeneous 
Neumann boundary conditions. The corresponding variational problem on V = 
H 1 ^) reads: 

Find u E V : a{v,u)=L(v) Vu E V, (2) 

where 

a(v,u) = / Vu • Vu + vudx, (3) 
Jn 

L(v) = [ vfdx. (4) 
Jn 

To assemble and solve a linear system AU = b for the degrees of freedom U G l w 
of a finite clement approximation Uh — ^2d—i Ui4>i E Vh C V, where the set of basis 
functions {4>%}f =1 spans V/,, one may simply define the bilinear (3) and linear (4) 
forms, and then call the two functions assemble and solve in DOLFIN. This is 
illustrated in Table I where we list a complete program for solving the reaction- 
diffusion problem (1) using piecewise linear elements. 

The example given in Table I illustrates the use of DOLFIN for solving a partic- 
ularly simple equation, but assembling and solving linear systems remain the two 
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key steps in the solution of more complex problems. We return to this in Section 5. 
2.1 Automated code generation 

DOLFIN may assemble a variational form of any rank 1 from a large class of vari- 
ational forms and it does so efficiently by automated code generation. Following a 
traditional paradigm, it is difficult to build automated systems that are at the same 
time general and efficient. Through automated code generation, one may build a 
system which is both general and efficient. 

DOLFIN relies on a form compiler to automatically generate code for the inner- 
most loop of the assembly algorithm from a high-level mathematical description of 
a finite element variational form, as discussed in Kirby and Logg [2006] and Olgaard 
et al. [2008]. As demonstrated in Kirby and Logg [2006], computer code can be 
generated which outperforms the usual hand-written code for a class of problems 
by using representations which can not reasonably be implemented by hand. Fur- 
thermore, automated optimization strategies can be employed [Kirby et al. 2005; 
Kirby et al. 2006; Kirby and Logg 2007; Olgaard and Wells 2010] and different 
representations can be used, with the most efficient representation depending on 
the nature of the differential equation [Olgaard and Wells 2010]. Recently, similar 
results have been demonstrated in SyFi/SFC [Alnaes and Mardal 2010]. 

Code generation adds an extra layer of complexity to a software system. For this 
reason, it is essential to isolate the parts of a program for which code must be gener- 
ated. The remaining parts may be implemented as reusable library components in a 
general purpose language. Such library components include data structures and al- 
gorithms for linear algebra (matrices, vectors and linear/nonlinear solvers), compu- 
tational meshes, representation of functions, input/output and plotting. However, 
the assembly of a linear system from a given finite element variational formulation 
must be implemented differently for each particular formulation and for each par- 
ticular choice of finite clement function space(s). In particular, the innermost loop 
of the assembly algorithm varies for each particular problem. DOLFIN follows a 
strategy of rc-usablc components at higher levels, but relies on a form compiler to 
generate the code for the innermost loop from a user-defined high-level description 
of the finite element variational form. 

DOLFIN partitions the user input into two subsets: data that may only be 
handled efficiently by special purpose code, and data that can be efficiently handled 
by general purpose library components. For a typical finite element application, 
the first set of data may consist of a finite element variational problem and the 
finite element(s) used to define it. The second set of data consists of the mesh 
and possibly other parameters that define the problem. The first set of data is 
given to a form compiler that generates special purpose code. That special purpose 
code may then use the second set of data as input to compute the solution. If the 
form compiler is implemented as a just-in-time (JIT) compiler, one may seamlessly 
integrate the code generation into a problem solving environment to automatically 
generate, compile and execute generated code at run-time on demand. We present 
this process schematically in Figure 2. 



x Rank refers here to the number of arguments to the form. Thus, a linear form has rank one, a 
bilinear form rank two, etc. 
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Fig. 2. An automated system (DOLFIN) using a JIT compiler to generate special purpose code 
for a subset of its input. For a typical finite element application, Input 1 consists of the variational 
problem and the finite elcment(s) used to define it, and Input 2 consists of the mesh. 



2.2 Compilation of variational forms 

Users of DOLFIN may use one of the two form compilers FFC or SyFi/SFC 2 
to generate problem-specific code. When writing a C++ application based on 
DOLFIN, the form compiler must be called explicitly from the command-line prior 
to compile-time. The form compiler generates CH — h code which may be included 
in a user program. The generated code defines a number of classes that may be 
instantiated by the user and passed to the DOLFIN C++ library. In particular, the 
user may instantiate form objects which correspond to the variational forms given 
to the form compiler and which may be passed as input arguments to the assembly 
function in DOLFIN. When using DOLFIN from Python, DOLFIN automatically 
handles the communication with the form compiler, the compilation (and caching) 
of the generated code and the instantiation of the generated form classes at run-time 
(JIT compilation). 

3. DESIGN CONSIDERATIONS 

The development of DOLFIN has been driven by two keys factors. The first is striv- 
ing for technical innovation. Examples of this include the use of a form compiler to 
generate code and new data structures for efficient representation of computational 
meshes [Logg 2009]. A second driving force is provided by the needs of applications; 
diverse and challenging applications have demanded and resulted in generic solu- 
tions for broad classes of problems. Often the canonical examples have not exposed 
limitations in the technology, particularly with respect to how the time required for 
code generation scales with the complexity of the considered equation [Olgaard and 
Wells 2010]. These have only become evident and then addressed when attempting 
to solve challenging problems at the limits of current technology. It is our expe- 
rience that both these components are necessary to drive advances and promote 
innovation. We comment in this section on some generic design considerations that 
have been important in the development of DOLFIN. 

3.1 Languages and language features 

DOLFIN is written primarily in CH — h with interfaces provided both in the form 
of a C++ class library and a Python module. The bulk of the Python interface 



2 DOLFIN may be used in conjunction with any form compiler conforming to the UFC interface. 
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is generated automatically using SWIG [SWIG 2009; Beazley 2003], with some 
extensions hand- written in Python 3 . The Python interface offers the performance 
of the underlying CH — h library with the ease of an intuitive scripting language. 
Performance critical operations are developed in C++, and users can develop solvers 
based on DOLFIN using either the C++ or Python interface. 

A number of C++ libraries for finite element analysis make extensive use of 
templates. Templated classes afford considerable flexibility and can be particularly 
useful in combining high-level abstractions and code re-use with performance as 
they avoid the cost inherent in virtual function calls in CH — h However, the extensive 
use of templates can obfuscate code, it increases compilation times and compiler 
generated error messages are usually expansive and difficult to interpret. We use 
templates in DOLFIN where performance demands it, and where it enables reuse of 
code. However, a number of key operations in a finite element library which require 
a function call involve a non-trivial number of operations within the function, and 
in these cases we make use of traditional C++ polymorphism. This enhances 
readability and simplifies debugging compared to template-based solutions, while 
not affecting run-time performance since the extra cost of a virtual function call 
is negligible compared to, for example, computing an element matrix or inserting 
the entries of an element matrix into a global sparse matrix. At the highest levels 
of abstraction, users are exposed to very few templated classes and objects, which 
simplifies the syntax of user-developed solvers. The limited use of templates at the 
user level also simplifies the automated generation of the DOLFIN Python interface. 

In mirroring mathematical concepts in the library design, sharing of data between 
objects has proved important. For example, objects representing functions may 
share a common object representing a function space, and different function spaces 
may share a common object representing a mesh. We have dealt with this issue 
through the use of shared pointers, and in particular boost : : sharebLptr from 
Boost. In managing data sharing, this solution has reduced the complexity of 
classes and improved the robustness of the library. While we make use of shared 
pointers, they are generally transparent to the user and need not be used in the 
high-level interface, thereby not burdening a user with the more complicated syntax. 

3.2 Interfaces 

Many scientific libraries perform a limited number of specialized operations which 
permits exposing users to a minimal, high-level interface. DOLFIN provides such 
a high-level interface for solving partial differential equations, which in many cases 
allows non-trivial problems to be solved with less than 20 lines of code (as we will 
demonstrate in Section 5). At the same time, it is recognized that methods for 
solving partial differential equations are diverse and evolving. Therefore, DOLFIN 
provides interfaces of varying complexity levels. For some problems, the minimal 
high-level interface will suffice, whereas other problems may be solved using a mix- 
ture of high- and low-level interfaces. In particular, users may often rely on the 
DOLFIN Function class to store and hide the degrees of freedom of a finite el- 



3 Thesc extensions deal primarily with JIT compilation, i.e., code generation, assembly and wrap- 
ping, of objects before sending them through the SWIG-generated Python interface to the under- 
lying C++ library. 
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ement function. Nevertheless, the degrees of freedom of a function may still be 
manipulated directly if desired. 

The high-level interface of DOLFIN is based on a small number of classes rep- 
resenting common mathematical abstractions. These include the classes Matrix, 
Vector, Mesh, FunctionSpace, Function and VariationalProblem. In addition 
to these classes, DOLFIN provides a small number of free functions, including 
assemble, solve and plot. We discuss these classes and functions in more detail 
in Section 4. 

DOLFIN relies on external libraries for a number of important tasks, including 
the solution of linear systems. In cases where functionality provided by external 
libraries must be exposed to the user, simplified wrappers are provided. This way, 
DOLFIN preserves a consistent user interface, while allowing different external 
libraries which perform similar tasks to be seamlessly interchanged. It also permits 
DOLFIN to set sensible default options for libraries with complex interfaces that 
require a large number of parameters to be set. This is most evident in the use of 
libraries for linear algebra. While the simplified wrappers defined by DOLFIN are 
usually sufficient, access is permitted to the underlying wrapped objects so that 
advanced users may operate directly on those objects when necessary. 

4. DESIGN AND IMPLEMENTATION 

Like many other finite element libraries, DOLFIN is designed as a collection of 
classes partitioned into components/libraries of related classes. However, while 
these classes are typically implemented as part of the library, see, e.g., Bangerth 
et al. [2007], DOLFIN relies on automated code generation and external libraries 
for the implementation of a large part of the functionality. Figure 3 shows a UML 
diagram of the central components and classes of DOLFIN. These include the lin- 
ear algebra classes, mesh classes, finite element classes and function classes. As 
already touched upon above, the linear algebra classes consist mostly of wrap- 
per classes for external libraries. The finite element classes Form, FiniteElement 
and Dof Map are also wrapper classes but for generated code, whereas the classes 
Assembler, VariationalProblem together with the mesh and function classes are 
implemented as regular C++ classes (with Python wrappers) as part of DOLFIN. 
In the following, we address these key components of DOLFIN, in order of increas- 
ing abstraction. In addition to the components depicted in Figure 3, DOLFIN 
includes a number of additional components for input /output, logging, plotting 
and the solution of ordinary differential equations. 

4.1 Linear algebra 

DOLFIN allows the transparent use of various specialized linear algebra libraries. 
This includes the use of data structures for sparse and dense matrices, precondi- 
tioners and iterative solvers, and direct linear solvers. This approach allows users 
to leverage the particular strengths of different libraries through a simple and uni- 
form interface. Currently supported linear algebra backends include PETSc [Balay 
et al. 2009], Trilinos/Epetra [Heroux et al. 2005], uBLAS 4 [Walter et al. 2009] and 



4 Krylov solvers and preconditioners for uBLAS are implemented as part of DOLFIN. 
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Fig. 3. UML diagram of the central components and classes of DOLFIN. 



MTL4 [Gottschling and Lumsdaine 2009]. Interfaces to the direct solvers UMF- 
PACK [Davis 2004] (sparse LU decomposition) and CHOLMOD [Chen et al. 2008] 
(sparse Cholesky decomposition) are also provided. 

The implementation of the DOLFIN linear algebra interface is based on C++ 
polymorphism. A common abstract base class GenericMatrix defines a minimal 
matrix interface suitable for finite element assembly, and a subclass of Generic- 
Matrix implements the functionality for each backend by suitably wrapping native 
data structures of its respective backend. Similarly, a common abstract base class 
GenericVector defines a minimal interface for vectors with subclasses for all back- 
ends. The two interface classes GenericMatrix and GenericVector are themselves 
subclasses of a common base class GenericTensor. This enables DOLFIN to im- 
plement a common assembly algorithm for all matrices, vectors and scalars (or any 
other rank tensor) for all linear algebra backends. Compared to a template-based 
solution, polymorphism may incur overhead associated with the cost of resolving 
virtual function calls. However, the most performance-critical function call to the 
linear algebra backend during assembly is typically insertion of a local element ma- 
trix into a global sparse matrix. This operation usually involves a considerable 
amount of computation/memory access, hence the extra cost of the virtual func- 
tion call in this case may be neglected. For cases in which the overhead of a virtual 
function call is not negligible, operating directly on the underlying object avoids 
this overhead. 

4.2 Meshes 

The DOLFIN Mesh class is based on a simple abstraction that allows dimension- 
independence, both in the implementation of the DOLFIN mesh library and in user 
code. In particular, the DOLFIN assembly algorithm is common for all simplex 
meshes in one, two and three space dimensions. We provide here an overview of 
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Table II. DOLFIN mesh abstractions and corresponding classes. Users may refer to a mesh entity 
either by a topological dimension and index or as a named mesh entity such as a vertex with a 
specific index. 

C++ code 

Mesh mesh ( "mesh. xml") ; 

for (Celllterator cell (mesh); Icell.endO; ++cell) 

for (Vertexlterator vertex(*cell) ; ! vertex. end() ; ++vertex) 
cout << vertex->dim() << " " << vertex->lndex() << endl; 



Python code 

mesh = MeshO'mesh.xml") 

for cell in cells (mesh): 

for vertex in vertices (cell) : 

print vertex. dim() , vertex. indexO 



Table III. Basic use of DOLFIN mesh iterators for iterating over all vertices of all cells of a mesh 
in C++ (top) and Python (bottom). 

the DOLFIN mesh implementation and refer to Logg [2009] for details. While only 
simplices are currently supported, the design paradigm extends to non-simplicial 
meshes. 

A DOLFIN mesh consists of a collection of mesh entities that define the topology 
of the mesh, together with a geometric mapping embedding the mesh entities in 
W 1 . A mesh entity is a pair (d,i), where d is the topological dimension of the 
mesh entity and i is a unique index of the mesh entity. A similar approach may 
be found in Knepley and Karpeev [2009]. Mesh entities are numbered within each 
topological dimension from to rid — 1, where rid is the number of mesh entities of 
topological dimension d. For convenience, mesh entities of topological dimension 
are referred to as vertices, entities of dimension 1 edges, entities of dimension 2 
faces, entities of codimension 1 facets and entities of codimension cells. These 
concepts are summarized in Table II. 

Algorithms operating on a mesh can often be expressed in terms of iterators [Berti 
2002; 2006]. The mesh library provides the general iterator MeshEntitylterator in 
addition to the specialized mesh iterators Vertexlterator, Edgelterator, Face- 
Iterator, Facetlterator and Celllterator. We illustrate the use of iterators in 
Table III. 
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The DOLFIN mesh library also introduces the concept of a mesh function. A 
mesh function, and its corresponding implementation MeshFunction, is a discrete 
function on the set of mesh entities of a specific dimension. It is only defined on a 
set of mesh entities which is in contrast to functions represented by the DOLFIN 
Function class which take a value at each point in the domain covered by the mesh. 
The class MeshFunction is templated over the value type which allows users, for 
example, to create a boolean- valued mesh function over the cells of a mesh to 
indicate regions for mesh refinement, an integer-valued mesh function on vertices 
to indicate a mapping from local to global vertex numbers or a float- valued mesh 
function on cells to indicate material data. 

The simple object-oriented interface of the DOLFIN mesh library is combined 
with efficient storage of the underlying mesh data structures. Objects like vertices, 
edges and faces are never stored. Instead, DOLFIN stores all mesh data in plain 
C/C++ arrays and provides views of the underlying data in the form of the class 
MeshEntity and its subclasses Vertex, Edge, Face, Facet and Cell, together with 
their corresponding iterator classes. An earlier version of the DOLFIN mesh li- 
brary used a full object-oriented model also for storage, but the simple array-based 
approach has reduced storage requirements and improved the speed of accessing 
mesh data by orders of magnitude [Logg 2009]. In its initial state, the DOLFIN 
Mesh class only stores vertex coordinates, using a single array of double values, 
and cell-vertex connectivity, using a compressed row-like data structure consisting 
of two arrays of unsigned int values. Any other connectivity, such as, vertex- 
vertex, edge-cell or cell-facet connectivity, is automatically generated and stored 
when required. Thus, if a user solves a partial differential equation using piecewise 
linear elements on a tctrahedral mesh, only cell-vertex connectivity is required and 
so edges and faces are not generated. However, if quadratic elements are used, 
edges are automatically generated and cubic elements will lead to a generation of 
faces as well as edges. 

In addition to efficient representation of mesh data, the DOLFIN mesh library 
implements a number of algorithms which operate on meshes, including adaptive 
mesh refinement (using a Rivara [1991]-type method), mesh coarsening and mesh 
smoothing. DOLFIN does not provide support for mesh generation, except for 
a number of simple shapes like squares, boxes and spheres. The following code 
illustrates adaptive mesh refinement in DOLFIN: 

C+ + code 

MeshFunction<bool> cell_markers(mesh, mesh. topology () .dimO) ; 

for (Celllterator cell (mesh); ! cell . end() ; ++cell) 
{ 

if (...) 

cell_markers [*cell] = true; 
else 

cell_markers [*cell] = false; 

} 

mesh.ref ine(cell_markers) ; 
mesh. smooth () ; 
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4.3 Finite elements 

DOLFIN supports a wide range of finite elements. At present, the following ele- 
ments are supported: 

(1) i/ 1 -conforming finite elements: 

(a) CG q , arbitrary degree continuous Lagrange elements. 

(2) iJ(div)-conforming finite elements: 

(a) RTq, arbitrary degree Raviart-Thomas elements [Raviart and Thomas 
1977]; 

(b) BDM g , arbitrary degree Brezzi-Douglas-Marini elements [Brezzi et al. 
1985]; and 

(c) BDFM 9 , arbitrary degree Brezzi-Douglas-Fortin-Marini elements [Brezzi 
et al. 1987]. 

(3) if (curl)-conforming finite elements: 

(a) NED g , arbitrary degree Nedelec elements (first kind) [Nedelec 1980]. 

(4) L 2 -conforming finite elements: 

(a) DG g , arbitrary degree discontinuous Lagrange elements; and 

(b) CRi, first degree Crouzeix-Raviart 5 elements [Crouzeix and Raviart 1973]. 

Arbitrary combinations of the above elements may be used to define mixed elements. 
Thus, one may for example define a Taylor-Hood element by combining a vector- 
valued P2 element with a scalar Pi clement. Arbitrary nesting is supported, thus 
allowing a mixed Taylor-Hood element to be used as a building block in a coupled 
problem which involves more than just the velocity and pressure fields. In Section 5, 
we demonstrate the use of mixed elements for the Poisson equation. Presently, 
DOLFIN only supports elements defined on simplices. This is not a technical 
limitation in the library design, but rather a reflection of current user demand. 

DOLFIN relies on a form compiler such as FFC for the implementation of finite 
elements. FFC in turn relies on FIAT for tabulation of finite element basis functions 
on a reference element. In particular, for any given element family and degree q from 
the list of supported elements, FFC generates C++ code conforming to a common 
interface specification for finite elements which is part of the UFC interface. Thus, 
DOLFIN does not include a library of finite elements, but relies on automated code 
generation, either prior to compilc-time or at run-time, for the implementation of 
finite elements. The generated code may be used for efficient run-time evaluation 
of finite element basis functions, derivatives of basis functions and evaluation of 
degrees of freedom (applying the functionals to any given function). However, 
these functions are rarely accessed by users as a user is not usually exposed to the 
details of a finite element beyond its declaration, and since DOLFIN automates 
the assembly of variational forms based on code generation for evaluation of the 
element matrix. Detailed aspects of automated finite element code generation can 
be found in Olgaard et al. [2008] for discontinuous elements and in Rognes et al. 
[2009] for ff(div) and ff(curl) elements. 



5 Crouzeix-Raviart elements are sometimes referred to as C°-noraconforming. 
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4.4 Function spaces 

The concept of a function space plays a central role in the mathematical formulation 
of finite element methods for partial differential equations. DOLFIN mirrors this 
concept in the class FunctionSpace. This class defines a finite dimensional function 
space in terms of a Mesh, a FiniteElement and a DofMap (degree of freedom map): 

C+ + code 

class FunctionSpace 
{ 

public: 
private : 

boost :: shared_ptr<const Mesh> _mesh; 

boost :: shared_ptr<const FiniteElement> .element ; 

boost :: shared_ptr<const DofMap> _dofmap; 

}; 



The mesh defines the domain, the finite element defines the local basis on each 
cell and the degree of freedom map defines how local function spaces are patched 
together to form the global function space. 

For some problems, finite element spaces are not appropriate and a "quadrature 
function space" can be used. In such a "space", functions can be evaluated at 
discrete points (quadrature points) but not elsewhere, and derivatives cannot be 
computed. This concept is discussed in 01gaard et al. [2008] and 01gaard and Wells 
[2010]. 

Incorporating the mathematical concept of function spaces in the library design 
provides a powerful abstraction, especially for sharing data in a transparent and 
simple fashion. In particular, several functions may share the same function space 
and thus the same mesh, finite element and degree of freedom mapping. 

4.5 Functions 

Functions on a finite element function space are implemented in DOLFIN in the 
form of the Function class. A Function is expressed as a linear combination of 
basis functions on a discrete finite element or quadrature space. The expansion 
coefficients (degrees of freedom) of the Function are stored as a (Generic) Vector: 

C+ + code 

class Function 
{ 

public: 
private : 

boost :: shared_ptr<const FunctionSpace> _function_space; 
boost: : shared_ptr<GenericVector> _vector; 

}; 



A Function may be evaluated at arbitrary points on a finite element mesh, used 
as a coefficient in a variational form, saved to file for later visualization or plotted 
directly from within DOLFIN. The DOLFIN Function class is particularly powerful 
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for supplying and exchanging data between different models in coupled problems, 
as will be demonstrated in Section 5. 

Evaluation of Functions at arbitrary points is handled efficiently using the GNU 
Triangulated Surface Library [GTS 2009]. With the help of GTS, DOLFIN locates 
which cell of the Mesh of the FunctionSpace contains the given point. The function 
value may then be computed by evaluating the finite element basis functions at 
the given point (using the FiniteElement of the FunctionSpace) and multiplying 
with the appropriate coefficients in the Vector (determined using the DofMap of 
the FunctionSpace). 

4.6 Expressions 

Many times, it is appropriate to express a coefficient in a variational problem by 
an expression or an algorithm for evaluating the coefficient at a given point, rather 
than expressing it as a linear combination of basis functions (as in the Function 
class). Such coefficients may be conveniently implemented using the Expression 
class. 

An Expression is defined by a user through overloading the Expression: :eval 
function. This functor construct provides a powerful mechanism for defining com- 
plex coefficients. In particular, the functor construct allows a user to attach data to 
an Expression. A user may, for example, read data from a file in the constructor of 
an Expression subclass which is then later accessed in the eval callback function. 

For the definition of functions given by simple expressions, like f(x) = sin(x) 
or g(x,y) = sin(x) cos(y), the DOLFIN Python interface provides simple and au- 
tomated JIT compilation of expressions. While the Python interface does allow 
a user to overload the eval function from Python 6 , this may be inefficient as the 
call to eval involves a callback from C++ to a Python function and this may be 
called repeatedly during assembly (once or more on each cell). However, JIT com- 
pilation avoids this by automatically generating, compiling, wrapping and linking 
C++ subclasses of the Expression class. 

An Expression may be evaluated at arbitrary points on a finite element mesh, 
used as a coefficient in a variational form, projected or interpolated into a finite 
element function space or plotted directly from within DOLFIN. Table IV illustrates 
use of the DOLFIN Function and Expression classes in Python. 

4.7 Variational forms 

DOLFIN allows general variational forms to be expressed in a form language that 
mimics mathematical notation. For example, consider the bilinear form of the 
standard Stokes variational problem. This may be conveniently expressed in the 
form language as illustrated in Table V. The form language allows the expression 
of general multilinear forms of arity p on the product space V } ] x Vfi x • ■ • x V£ of 
a sequence {V^}j =1 of finite element spaces on a domain del", 

a : V£ x Vl x ■ ■ ■ x V h p -> R. (5) 

Such forms are fundamental building blocks in linear and nonlinear finite element 
analysis. In particular, linear forms (p = 1) and bilinear forms (p = 2) are central 



6 SWIG supports cross-language polymorphism using the director feature. 
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Python code 

# Create mesh 

mesh = UnitSquare (32 , 32) 

# Define an expression 

f = Expression(("sin(x[0])", "cos (x [1] ) ") ) 

# Project expression to a finite element space 
V = VectorFunctionSpaceCmesh, "CG", 2) 

g = project(f , V) 

# Evaluate expression and function 
print f(0.1, 0.2) 

print g(0.1, 0.2) 

# Plot expression and function 
plot(f, mesh=mesh) 

plot(g) 



Table IV. Defining, projecting, evaluating and plotting expressions and functions using the 
DOLFIN Python interface. 



Mathematical notation 



a(v, u) = /q grad v ■ grad u — div vp -\- q div u dx 



Code 



a = (inner (grad(v) , grad(u)) - div(v)*p + q*div(u))*dx 



Table V. Expressing the bilinear form for the Stokes equations in DOLFIN. 

to the finite clement discretization of partial differential equations. Forms of higher 
arity are also supported as they may sometimes be of interest, see Kirby and Logg 
[2006]. 

DOLFIN relies on the Unified Form Language (UFL) [Alnses and Logg 2009; 
Alnses 2009] for the expression of variational forms. The form language allows the 
expression of a wide range of finite element variational forms in a language close 
to mathematical notation. UFL also supports functional differentiation of general 
nonlinear forms. Forms can involve integrals over cells, interior facets and exterior 
facets. Line and surface integrals which do not coincide with cell facets are not yet 
supported, although extensions in this direction for modeling crack propagation 
have been made [Nikbakht and Wells 2009]. For details on the form language, we 
refer to the UFL user manual [Alnaes and Logg 2009]. 

A user of the DOLFIN C++ interface will typically define a set of forms in one or 
more form files and call FFC on the command-line. The generated code may then be 
included in the user's C++ program. As an illustration, consider again the bilinear 
form of the Stokes problem as expressed in Table V. This may be entered together 
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with the corresponding linear form L = v*f*dx in a text file named Stokes. ufl 
which may then be compiled with FFC: 



ffc -1 dolfin Stokes. ufl 



This will generate a C++ header file Stokes .h which a user may include in a C++ 
program to instantiate the pair of forms: 

C+ + code 

#include <dolfin.h> 
#include "Stokes. h" 

int mainO 
{ 

Stokes: : FunctionSpace V(mesh); 
Stokes: : BilinearForm a(V, V); 
Stokes: : LinearForm L(V) ; 

} 



When used from Python, form compilation is handled automatically by DOLFIN. If 
a form is encountered during the execution of a program, the necessary CH — h code 
is automatically generated and compiled. The generated object code is cached so 
that code is generated and compiled only when necessary. Thus, if a user solves 
the Stokes problem twice, code is only generated the first time, as the JIT compiler 
will recognize the Stokes form on subsequent runs. 

4.8 Finite element assembly 

Given a variational form, the DOLFIN assemble function assembles the corre- 
sponding global tensor. In particular, a matrix is assembled from a bilinear form, 
a vector is assembled from a linear form, and a scalar value is assembled from a 
rank zero form (a functional). While DOLFIN does not provide data structures 
for sparse tensors of rank greater than two, the abstract GenericTensor interface, 
which was introduced in Section 4.1, permits users to supply data structures for 
arbitrary rank tensors. 

To discretize the multilinear form (5), we may introduce a basis {4> J k }^l 1 for each 
function space V£, j = 1, 2, . . . , p, and define the global tensor 

A i = a(4,^ 2 , ••-,<), (6) 

where i — (ii, 12, ■ ■ ■ , i p ) is a multi-index. If the multilinear form is defined as 
an integral over O = LiKeT h K 7 the tensor A may be computed by assembling the 
contributions from all elements, 

A i = a(4> 1 il ,fc,...,ti p )= <^ 2 , ( 7 ) 

KET 

where a K denotes the contribution from element K. We further let {<^^Yk=\ 
denote the local finite element basis for V£ on K and define the element tensor A K 
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(the "element stiffness matrix") by 

A? = a^«' 1 , (8) 

The assembly of the global tensor A thus reduces to the computation of the element 
tensor A K on each element K and the insertion of the entries of A K into the 
global tensor A. In addition to contributions from all cells, DOLFIN also assembles 
contributions from all exterior facets (facets on the boundary) and all interior facets 
if required. 

The key to the generality and efficiency of the DOLFIN assembly algorithm 
lies in the automated generation of code for the evaluation of the element tensor. 
DOLFIN relies on generated code both for the evaluation of the element tensor 
and the mapping of degrees of freedom. Thus, the assembly algorithm may call the 
generated code on each cell of the mesh, first to compute the element tensor and then 
again to compute the local-to-global mapping by which the entries of the element 
tensor may be inserted into the global tensor. The complexity inherent in non- 
trivial forms, such as those which involve mixed element spaces, vector elements and 
discontinuous Galerkin methods, is not exposed in the form abstraction. DOLFIN 
is unaware of how the element matrix is represented or how forms are integrated. 
It simply provides coefficient and mesh data to the generated code and assembles 
the computed results. The algorithm for computing the element tensor is instead 
determined by the form compiler. Various algorithms are possible, including both 
quadrature and a special tensor representation, and the most efficient algorithm 
can depend heavily on the nature of the form [Kir by and Logg 2006; Olgaard and 
Wells 2010]. 

To assemble a Form a, a user may simply call the function assemble which 
computes and returns the corresponding tensor. Thus, a bilinear form may be 
assembled by 



C++ code 


Matrix A; 




assemble (A, a) ; 




C++ and 






Python code 


A = assemble (a) 



in Python. Several optional parameters may be given to specify either assembly 
over specific subdomains of the mesh or reuse of tensors. 

4.9 Boundary conditions 

Natural boundary conditions are enforced weakly as part of a variational problem 
and are typically of Neumann or Robin type, but may also be of Dirichlet type 
as will be demonstrated Section 5. Essential boundary conditions are typically of 
Dirichlet type and are enforced strongly at the linear algebra level. DOLFIN also 
supports the specification of periodic boundary conditions. We describe here the 
definition and application of strong Dirichlet boundary conditions. 

ACM Transactions on Mathematical Software, Vol. V, No. N, Month 20YY. 



18 • A. Logg and G. N. Wells 



We define a Dirichlet boundary condition in terms of a function space V, a 
function g and a subset of the boundary T C <9f2, 

u(x) = g(x) Vx e r. (9) 

The corresponding definition in the DOLFIN Python interface reads 

Python code 

be = DirichletBC(V, g, gamma) 



where V is a FunctionSpace, g is a Function or Expression, and gamma is a 
SubDomain. Alternatively, the boundary may be defined in terms of a MeshFunction 
marking a portion of the facets on the mesh boundary. The function space V defines 
the space to which the boundary condition will be applied. This is useful when 
applying a Dirichlet boundary condition to particular components of a mixed or 
vector-valued problem. 

Once a boundary condition has been defined, it can be applied in one of two 
ways. The simplest is to act upon the assembled global system: 

Python code 

be. apply (A, b) 



For each degree of freedom to be constrained, this call will zero the corresponding 
row in the matrix, set the diagonal entry to one and put the Dirichlet value at 
the corresponding position in the right-hand side vector. An optional argument 
can be provided for updating the boundary conditions inside a Newton iteration. 
Alternatively, the boundary condition may be supplied directly to the assembler 
which will then apply the boundary condition by modifying the element matrices 
in a manner that preserves any symmetry of the global matrix: 

Python code 

A, b = assemble_system(a, L, be) 



4.10 Variational problems 

At the highest level of abstraction, objects may be created that represent variational 
problems of the canonical form (2). Such a variational problem may be defined and 
solved by 

Python code 

problem = VariatlonalProblem(a, L) 
u = problem. solve () 



A constraint on the trial space in the form of one or more Dirichlet conditions may 
be supplied as additional arguments. Other parameters include the specification 
of the linear solver and preconditioner (when appropriate) and whether or not the 
variational problem is linear. In the case of a nonlinear variational problem where 
one seeks to satisfy 



F(v) = V v 6 V, 
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the bilinear form a is interpreted as the Gateaux derivative of a nonlinear form 
L = F. 

4.11 File I/O and visualization 

DOLFIN provides input/output for objects of all its central container classes, in- 
cluding Vector, Matrix, Mesh and MeshFunction. Objects are stored to file in 
XML format. For example, a Mesh may be loaded from and stored back to file by 

C+ + code 

File f ile("mesh.xml") ; 
Mesh mesh; 

file >> mesh; 
file << mesh; 



Mesh data may be converted to the native DOLFIN XML format from Gmsh, 
Medit, Diffpack, ABAQUS, Exodus II and StarCD formats using the conversion 
utility dolf in-convert. 

Solution data may be exported in a number of formats, including the VTK XML 
format which is useful for visualizing a Function in VTK-based tools, such as 
ParaView. DOLFIN also provides built-in plotting for Mesh, MeshFunction and 
Function using Viper [Skavhaug 2009] by 

C++ code 

plot (mesh) ; 

plot (mesh_f unction) ; 

plot (u) ; 



5. APPLICATIONS 

We present here a collection of examples to demonstrate the use of DOLFIN for solv- 
ing partial differential equations and related problems of interest. A more extensive 
range of examples are distributed with the DOLFIN source code. For a particularly 
complicated application to reservoir modeling, we refer to Wells et al. [2008]. Some 
issues of particular relevance to solid mechanics problems, such as plasticity, are 
discussed in Olgaard et al. [2008]. All examples correspond to DOLFIN 0.9.5. 

5.1 Evaluating functionals 

We begin with the simplest form that we can evaluate, a functional. In the absence 
of viscous stresses, the lift acting on a body can be computed by integrating the 
pressure multiplied by a suitable component of the unit vector normal to the surface 
of interest. The definition of this functional is shown in Table VI. From this 
definition, C++ code may be generated using a form compiler and then used to 
compute the lift generated by a computed pressure field. 

Another common application of functionals is the evaluation of various norms or 
evaluating the error of a computed solution when the exact solution is known. For 
example, one may define the squared L 2 and Hq norms as v*v*dx and dot (grad(v) , 
grad(v))*dx respectively. Alternatively, one may use the built-in DOLFIN func- 
tions norm and errornorm to evaluate norms and errors: 

ACM Transactions on Mathematical Software, Vol. V, No. N, Month 20YY. 



20 • A. Logg and G. N. Wells 



Form compiler code 



element = FiniteElementO'Lagrange" , triangle, 1) 



p = Function(element) 
n = triangle. n 
M = p*n[l]*ds 



Table VI. Definition of the functional for computing the lift. 



Python code 



print norm(v, "L2") 

print norm(v, "HI") 

print norm(v, "H10") 

print norm(v, "Hdiv") 

print errornorm(u_h, u, "L2") 

print errornorm(u_h, u, "HI") 

print errornorm(u_h, u, "H10") 

print errornorm(u_h, u, "Hdiv") 



5.2 Solving linear partial differential equations 

To illustrate the use of DOLFIN for solving simple linear partial differential equa- 
tions, we consider Poisson's equation — Au = / discrctized using three different 
methods: an H 1 -conforming primal approach using standard continuous Lagrange 
basis functions; a mixed method using H(div) x L 2 -conforming elements; and a 
discontinuous Galerkin method using L 2 -conforming Lagrange basis functions. 

5.2.1 H 1 -conforming discretization of Poisson's equation. For the standard H 1 - 
conforming approach, the bilinear and linear forms are given by 



V = FunctionSpace(mesh, "CG", 1) 
v = TestFunction(V) 
u = TrialFunction(V) 
f = Expression( . . . ) 

a = dot(grad(v), grad(u))*dx 
L = v*f*dx 



5.2.2 H(div)xL 2 - conforming discretization of Poisson's equation. For the mixed 
version of the Poisson problem, with u = on dil, the bilinear and linear forms 




(11) 



(12) 



and the forms may be specified in DOLFIN by 



Python code 
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read [Brezzi and Fortin 1991]: 



a(r, w; a, u) = / t • a — (V • r) u + w (V • a) Ax, (13) 
Jn 

L(t,w) = / wfdx, (14) 
Jn 

where t,<j £ V, w,u £ W and 

V = {t E H (div, fi) : r| K £ [P,] n (if) VA} , (15) 
lf-{t«eL 2 (fi) : ttflx G Pg_i (if) VA"} . (16) 

The corresponding implementation in DOLFIN for q = 2 reads: 

Python code 

V = FunctionSpace(mesh, "BDM", 2) 
W = FunctionSpace(mesh, "DG" , 1) 

mixed_space = V + W 

(tau, w) = TestFunctions(mixed_space) 
(sigma, u) = TrialFunctions(mixed_space) 

f = Expression( . . . ) 

a = (dot(tau, sigma) - div(tau)*u + w*div (sigma) ) *dx 
L = w*f*dx 

5.2.3 L 2 -conforming discretization of Poisson's equation. For a discontinuous 
interior penalty formulation of the Poisson problem, the bilinear and linear forms 
read: 

a(v,u)= / Vv-Vudx- [t>] • (Vu) ds - / (Vv) ■ [u] ds 
Jn\r° Jr° Jr° 

— / wn-Vuds— / \7v-unds+ / — fvj ■ [«] ds + / j-vuds (17) 
Jan J on Jt° h Jan " 

and 

L(v)= f vfdx, (18) 

where T° denotes all interior facets and v,u eV = {v e L 2 (Q) : v\ K e P q (K) \/K). 
Using ds to denote integration over exterior facets and dS to denote integration over 
interior facets, the corresponding implementation in DOLFIN for q = 1 reads as 
follows: 
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Python code 



V = FunctionSpace(mesh, "DG", 1) 

v = TestFunction(V) 
u = TrialFunction(V) 
f = Expression( . . . ) 

n = FacetNormal(mesh) 
h = CellSize (mesh) 

alpha =4.0 

a = dot (grad(v) , grad(u))*dx \ 

- dot (jump(v, n) , avg(grad(u) ) ) *dS \ 

- dot (avg(grad(v) ) , jump(u, n))*dS \ 

- v*dot (grad(u) , n)*ds - dot(grad(v), n)*u*ds \ 
+ alpha/h('+ , )*dot(jump(v, n) , jump(u, n))*dS \ 
+ (alpha/h)*v*u*ds 

L = v*f*dx 



5.3 Solving time-dependent partial differential equations 

Unsteady problems can be solved by defining a variational problem to be solved in 
each time step. We illustrate this by solving the convection-diffusion problem 



The velocity field b = b(x) may be a user-defined expression or an earlier computed 
solution. Multiplying (19) with a test function and discretizing in time using the 
Crank-Nicolson method, we obtain 



where k n — t n — t n -i is the time step and x — (x n + x n ~ 1 )/2. We may implement 
the problem (20) in DOLFIN by moving all terms involving to the right-hand 
side. Alternatively we may rely on the built-in operators lhs and rhs to extract 
the pair of bilinear and linear forms as illustrated in Table VII. In Table VIII we 
show the corresponding C++ program. 

5.4 Solving nonlinear partial differential equations 

Solution procedures for nonlinear differential equations are inherently more complex 
and diverse than those for linear equations. With this in mind, the design of 
DOLFIN allows users to build complex solution algorithms for nonlinear problems 
using the basic building blocks assemble and solve. However, a built-in Newton 
solver is also provided which suffices for many problems. We illustrate the solution 
of a nonlinear problem for the following nonlinear Poisson-like equation: 



u + b ■ Vu - V ■ (cVu) = /. 



(19) 




(20) 



V • (1 + u 2 ) V« = / in O, 



(21) 
(22) 



u = on dft. 
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Form compiler code 



scalar = FiniteElement ( "Lagrange " , triangle, 1) 
vector = VectorElement ( "Lagrange " , triangle, 2) 



v = TestFunction(scalar) # 

ul = TrialFunction(scalar) # 

uO = Function(scalar) # 

b = Function(vector) # 

f = Function(scalar) # 

c = 0.005 # 

k = 0.05 # 



test function 
solution at t_n 
solution at t_{n-l} 
convective velocity 
source term 
dif f usivity 
time step 



u = 0.5*(u0 + ul) 

F = v*(ul - u0)*dx + k*v*dot(b, grad(u) ) *dx + k*c*dot (grad(v) , grad(u))*dx 
a = lhs(F) 

L = rhs(F) + k*v*f*dx 



Table VII. Specification of the variational problem for the unsteady convection-diffusion equa- 
tion (19). 

Multiplying by a test function v e V = Hq(Q) and integrating over the domain f2, 
we obtain 



F(v;u)= / (l + u 2 ) Vv Vudx- / vf = 0, 
Jn Jn 



(23) 



where we note that F : V x V — > R is linear in its first argument and nonlinear 
in its second argument. To solve the nonlinear problem by Newton's method, we 
compute the Gateaux derivative D u F(v; u) and obtain 

dF(v; u + eSu) 



a(v, Su; it) = D u F(v; u)Su 



de 



e=0 (24) 

= / (f + u 2 ) Vv- \76udx+ / 2u5uVv ■ Vu dor. 
Jn Jn 

We note that F(-; u) : V — >• M is a linear form for every fixed u and that a(-, •; u) : 
V x V — > K is a bilinear form for every fixed u. A full solver for (21)-(22) in the case 
f(x,y) = xsiny is presented in Table IX. The form language UFL supports auto- 
matic differentiation, so many problems, including this one, can also be linearized 
automatically. 

6. CONCLUSIONS 

We have presented a problem solving environment that largely automates the fi- 
nite element approximation of solutions to differential equations. This is achieved 
by generating computer code for parts of the problem which are specific to the 
considered differential equation, and designing a generic library which reflects the 
mathematical structure of finite element variational problems. Using a high level 
of mathematical abstraction and automated code generation, the system can be 
designed for both readability and performance, allowing new models to be imple- 
mented rapidly and solved efficiently. 

Until recently, the focus has been on automating the assembly of linear systems 
arising from the finite element discretization of variational problems, in particular 
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C+ + code 



1/ Read mesh from file 
Mesh meshC'mesh.xml.gz") ; 

// Read velocity field from file 
Velocity: : FunctionSpace W(mesh) ; 
Function velocity(W, " velocity. xml. gz" ) ; 

// Read sub domain markers from file 

MeshFunction<unsigned int> sub_domains(mesh, "sub-domains. xml. gz") ; 
// Create function space 

ConvectionDif fusion: : FunctionSpace V(mesh) ; 

// Create source term and initial condition 
Constant f(0); 
Function u(V) ; 

// Set up variational forms 

ConvectionDif fusion: :BilinearForm a(V, V); 

a.b = velocity; 

ConvectionDif fusion: :LinearForm L(V); 
L.uO = u; L.b = velocity; L.f = f; 

// Set up boundary condition 
Constant g(l); 

DirichletBC bc(V, g, sub_domains, 1); 

// Linear system 
Matrix A; 
Vector b; 

// Assemble matrix and apply boundary conditions 
assemble (A, a); 
be. apply (A) ; 

// Parameters for time-stepping 

double T = 2.0; double k = 0.05; double t = k; 

// Output file 

File f ile( "temperature .pvd") ; 

// Time-stepping 
while (t < T) 
{ 

assemble (b, L) ; 
be. apply (b) ; 

solve(A, u.vectorQ, b, lu) ; 
file « u; 
t += k; 

} 



Tabic VIII. Implementation of the solver for the unsteady convection-diffusion equation 
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Python code 

from dolfin import * 

# Create mesh and define function space 
mesh = UnitSquare (32 , 32) 

V = FunctionSpace(mesh, "CG", 1) 

# Define boundary condition 

be = DirichletBC(V, Constant(O), DomainBoundary ( ) ) 

# Define source term and solution function 
f = Expression("x[0] *sin(x[l] ) ") 

u = Function(V) 

# Define variational problem 
v = TestFunction(V) 

du = TrialFunction(V) 

a = (1.0 + u*u)*dot(grad(v) , grad(du))*dx + \ 

2*u*du*dot (grad(v) , grad(u))*dx 
L = (1.0 + u*u)*dot(grad(v) , grad(u))*dx - v*f*dx 

# Solve nonlinear variational problem 

problem = VariationalProblem(a, L, be, nonlinear=True) 
problem . solve (u) 

# Plot solution and solution gradient 
plot (u) 

plot(grad(u)) 



Table IX. Implementation of a solver for the nonlinear Poisson problem (21)-(22). 

with regards to providing a general implementation independent of the variational 
problem, the mesh, the discretizing finite element space(s) and the linear algebra 
backend. More recently, efficient parallel computing has been added and automated 
error estimation/adaptivity is being developed. 
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